% This is the dynare code for Technical default model without covenants
% Authors: Emilio Bisetti, Kai Li, and Jun Yu
% Date: 2024, May 14
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


% Endogenous variables
var 
% 18 Endogenous variables
y   % output
ii   % Investment 
rl  % loan rate
rd  % risk free rate on deposit
rk  % return on capital
q   % price of purchasing capital
m   % Household's SDF
n   % Net worth
b   % debt
act_def % actual default prob
spread  % loan spread
re   % equity return
premiumH % external finance premium
c   % Consumption
lamda   % Marginal value of net worth
omegahat % default cutoff
h    % the debt-to-networth ratio
dk   % Growth rate of aggregate capital
a    % Exogenous TFP shock
dh   % household deposit
sigma % standard deviation of omega shock
u   % Household utility value
eu  % Household's expected utility value

%% 7 auxiliary variables
dn   % Growth rate of net worth
dc   % Growth rate of consumption
di   % Growth rate of investment
dy   % Growth rate of output
db   % Growth rate of debt
dq   % Growth rate of capital price
g    % Capital adjustment cost
;
  
% in total 25 endo variables  

% Exogenous variables
varexo ea;

% Parameters 

parameters 
ppsi    % Relative risk aversion
beta    % discount rate
gamma   % Intertempral elasticity of substitution
a1      % capital adjustment cost function parameter
xi      % capital adjustment cost function parameter
a2      % capital adjustment cost function parameter
alppha  % capital share
delta   % capital depriciation rate
rho_a   % Persistence of TFP shock
sigma_a % STD of TFP shock
eta     % Entrepreneur servival rate
vphi    % bank's verfication or mornitoring cost
chi     % fraction of redrawed networth
mu_ss  % mean of omega shock
sigma_ss % steady state of standard deviatio of omega
rho_s  % Persistence of omega1 shock
sigma_s % STD of sigma1 shock

%% The steady state value as initial value
yss    % output
iiss   % Investment
rlss   % loan rate
rdss   % risk free rate on deposit
rkss   % return on capital
qss    % price of purchasing capital
mss    % Household's SDF
nss    % Net worth
bss    % debt
uss    % Household utility value
euss   % Household's expected utility value
css    % Consumption
lamdass   % Marginal value of net worth
dkss   % Growth rate of aggregate capital
ass    % Exogenous TFP shock
dhss   % household deposit
omegahatss  % the default cutoff
hss    % the debt-to-networth ratio
ress   % return on equity
act_defss % acutal default probability
;        

load para_v1;

set_param_value('ppsi',ppsi)
set_param_value('beta',beta)
set_param_value('gamma',gamma)
set_param_value('a1',a1)
set_param_value('xi',xi)
set_param_value('a2',a2)
set_param_value('alppha',alppha)
set_param_value('delta',delta)
set_param_value('ass',log(Ass))
set_param_value('rho_a',rho_a)
set_param_value('rho_s',rho_s)
set_param_value('sigma_s',sigma_s)
set_param_value('sigma_a',sigma_a)
set_param_value('eta',eta)
set_param_value('vphi',zeta)
set_param_value('chi',chi)
%set_param_value('mu_ss',mu_sum)

%% steady state values
set_param_value('sigma_ss',log(sig_sum))
set_param_value('yss',log(Yss))
set_param_value('iiss',log(Iss))
set_param_value('rlss',log(RBss))
set_param_value('ress',log(REss))
set_param_value('rdss',log(RDss))
set_param_value('rkss',log(RKss))
set_param_value('qss',log(Qss))
set_param_value('mss',log(Mss))
set_param_value('nss',log(Nss))
set_param_value('bss',log(Bss))
set_param_value('css',log(Css))
set_param_value('lamdass',log(lamdass))
set_param_value('dkss',log(DK))
set_param_value('ass',log(Ass))
set_param_value('dhss',log(DHss))
set_param_value('hss',log(Hss))
set_param_value('uss',uss)
set_param_value('euss',euss)
set_param_value('act_defss',act_defss)
set_param_value('omegahatss',omegahatss)


model;
%% Household Side
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% dc is actually the actual growth rate of unnormalized C
% dy, di, dj is of the same form
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Auxiliary expressions.  These simplify the equations without adding
% additional variables.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
# z          = (log(omegahat) + exp(2*sigma(-1)) / 2) / exp(sigma(-1));
# zplus      = (log(omegahat(+1)) + exp(2*sigma) / 2) / exp(sigma);
# F          = normcdf(z);
# G          = normcdf(z - exp(sigma(-1)));
# GAMMA      = omegahat*(1-F)+G;
# Fp1        = normcdf(zplus);
# Gp1        = normcdf(zplus - exp(sigma));
# GAMMAp1    = omegahat(+1)*(1-Fp1)+Gp1;
# dGAMMA     = 1-F;
# dF         = normpdf(z)/(omegahat*exp(sigma(-1)));
# dG         = omegahat*dF;
# dFp1       = normpdf(zplus)/(omegahat(+1)*exp(sigma));
# dGAMMAp1   = 1-Fp1;
# dGp1       = omegahat(+1)*dFp1;

% 1 HHs preference 
exp(u*(1-1/ppsi)) = (1-beta)*exp((1-1/ppsi)*c) + beta*exp(eu*(1-1/ppsi));

% 2 Auxiliary eu
exp((1-gamma)*eu) = exp((1-gamma)*(u(+1) + dk));

% 3 SDF m
m = log(beta) + (-1/ppsi)*(c - c(-1) + dk(-1)) + (1/ppsi - gamma)*(u - eu(-1) + dk(-1));

% 4 HH's Euler equation
1/exp(rd) = exp(m(+1));

% 5 Return on capital Rk
exp(rk) = (alppha* exp(a) + exp(q)*(1-delta))/exp(q(-1));

%% Capital good producer
% 6 capital adjustment cost
g = ln(a1/(1-1/xi)*exp((1-1/xi)*ii)+a2);

% 7 Capital price
q = 1/xi*ii- ln(a1);

% 8 Production Function
exp(y) = exp(a);

% 9 Law of motion of capital
exp(dk) = (1-delta) + exp(g);

% 10 TFP shock
a - ass = rho_a*(a(-1)-ass) - sigma_a*ea;

%% Entrepreneurs problem
% 11 Balance sheet condition
exp(n) + exp(b) = exp(q+dk);


% 12 Debt
exp(b)=exp(n)*exp(h);

% 13 debt market clearing condition
exp(b) = exp(dh);

%% 14 Market Clearing condition
exp(y) = exp(c) + exp(ii)+ vphi*exp(n(-1)+rk)*(1+exp(h(-1)))*G;

% 15 FOC of H 
exp(m(+1)+rk(+1))*(1-GAMMAp1) + exp(lamda)*exp(m(+1)+rk(+1))*(GAMMAp1-vphi*Gp1) = exp(lamda);

% 16 FOC of omegahar
exp(lamda)=exp(m(+1)+rk(+1))*dGAMMAp1/(exp(m(+1)+rk(+1))*(dGAMMAp1-vphi*dGp1));

% 17 lender's break even condition
exp(h)=exp(m(+1)+rk(+1))*(GAMMAp1-vphi*Gp1)*(1+exp(h));

% 18 Law of motion of net worth
exp(n)*exp(dk(-1)) = exp(n(-1))*eta*(exp(rk)*(1+exp(h(-1))))*(1-GAMMA)+(1-eta*(1-F))*chi*exp(n(-1));

% 19 the actual default prob
act_def = F;

% 20 loan spread
spread = rl-rd;

% Auxiliary Growth rates 
% 21
dc = c - c(-1) + dk(-1);
% 22
dy = y - y(-1) + dk(-1);
% 23
di = ii - ii(-1) + dk(-1);
% 24
dn = n - n(-1) + dk(-1);
% 25
db = b - b(-1) + dk(-1);
% 26 
dq = q - q(-1);

% 27 return on equity
exp(re) = exp(rk)*(1+exp(h(-1)))*(1-GAMMA);

% 28 External finance premium        
premiumH = rk(+1) - rd; % definiton of primium, Rk/R

%29 sigma Shock 
sigma - sigma_ss = rho_s*(sigma(-1)-sigma_ss)+sigma_s*ea;


%30 loan rate equation
omegahat*exp(rk+q(-1)+dk(-1)-b(-1))=exp(rl(-1));

end;          
% in total 29 obs equations
steady_state_model;
%initval;
y =yss;  % output
ii = iiss;  % Investment
rl =rlss; % loan rate
rd =rdss; % risk free rate on deposit
re =ress; % equity return
rk =rkss; % return on capital
premiumH = rk-rd;
q  =qss; % price of purchasing capital
m  =mss; % Household's SDF
n =nss;  % Net worth
b =bss;  % debt
u =uss;  % Household utility value
eu =euss; % Household's expected utility value
c  =css; % Consumption
lamda =lamdass;  % Marginal value of net worth
dk =dkss;  % Growth rate of aggregate capital
a =ass; % Exogenous TFP shock
dh = dhss; % household deposit
omegahat = omegahatss; % the default cutoff
h =hss ; % the debt-to-networth ratio
g = iiss;
dc = dkss;
dy   = dkss;
di   = dkss;
dn   = dkss;
db = dkss;
dq = 0;
sigma = sigma_ss;
act_def=act_defss;
spread = rl-rd;

end;

shocks;
var ea; stderr 1;
end;


steady;
check;
model_diagnostics;

write_latex_dynamic_model;
% set_dynare_seed(1000);
% options_.relative_irf
options_.pruning = 1;
stoch_simul(order=2,periods=2000,drop=800,irf=30, pruning,replic=1000,simul_replic=1000); 

%read out simulations
simulated_series_raw=get_simul_replications(M_,options_);
save('simulated_series_raw.mat','simulated_series_raw');


